Thermodynamic functions of harmonic Coulomb crystals 
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Phonon frequency moments and thermodynamic functions (electrostatic and vibrational parts of 
the free energy, internal energy, and heat capacity) are calculated for bcc and fee Coulomb crystals 
in the harmonic approximation with a fractional accuracy < 10 -5 . Temperature dependence of 
thermodynamic functions is fitted by analytic formulas with an accuracy of a few parts in 10 5 . 
The static-lattice (Madelung) part of the free energy is calculated with an accuracy ~ 10~ 12 . The 
Madelung constant and frequency moments of hep crystals are also computed. 

<D . PACS numbers: 52.27.Gr, 52.25.Kn, 05.70.Ce, 97.20.Rp 
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I. INTRODUCTION 
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Coulomb crystals, introduced into theory by Wigncr have been studied by many authors. A thorough discussion 
was given, e.g., in Refs. (|||. Ewald technique has been used to calculate shear constants (6| and dispersion 
I of such crystals. The thermodynamics in the harmonic-lattice approximation has been analyzed, e.g., 
in Refs. |^|,||. Anharmonic corrections have been discussed in Refs. JlO)— [l2|] . Chabrier et al. jl3| have suggested 
' o ■ an approximate analytic model of the harmonic Coulomb crystal, which is widely used in astrophysics (e.g., Refs. 

P, ^4p^ ]). However, precise numerical calculations of the thermodynamic functions, valid at any temperature T, have 
, not been published. 

Here we report highly accurate calculations of phonon spectra and frequency moments of body-centered-cubic 
J>^, (bcc), face-centered-cubic (fee), and hexagonal-close-packed (hep) one-component Coulomb lattices in the harmonic 
approximation. We present also accurate calculations of thermodynamic functions for bcc and fee lattices at any values 
of the quantum parameter — T p /T, where T p = hLOp/ks is the ion plasma temperature and uj p — y 4TrniZ 2 e 2 jM 
is the ion plasma frequency (m, M, and Ze being the ion number density, mass, and charge, respectively). The 
' numerical results are given in the easy-to-use form of tables and fitting formulas. 



II. PHONON SPECTRUM AND ELECTROSTATIC ENERGY 



Consider a crystal of identical ions immersed in the uniform compensating background. The basic definitions are 
as follows (e.g., Ref. ||). Let us take an arbitrary ion as the origin of a Cartesian reference frame and specify the 
lattice basis li, I2, 13 generating direct lattice vectors I(ni, ri2, rc.3) — + + n^ls, where rii, n^, are arbitrary 
integers. The vectors g(«i, «2, "3) = ™igi + ™2g2 + ^3g3j where g^ ■ L, = 27T<%, form the reciprocal lattice. Consider 
also the primitive cell, the parallelepiped {i^li + ^2 + ^313} with < v\, u%, 1/3 < 1. Let 2V ce u be the number of ions 
in the primitive cell enumerated with an index k. The choice of the vectors 1^ is not unique, and one can describe 
a given lattice using different Y cc n. We will adopt the standard convention and choose the primitive cell with the 
lowest iV cc ii. The bcc and fee lattices are simple (the lowest Y co n = 1), whereas for the hep lattice one has the lowest 
C^h; AUi = 2. 

Along with the primitive cell one usually considers the Wigner-Seitz (WS) cell, which is a polyhedron with faces 
• '-j i crossing the lattice vectors at their midpoints at the right angle. The volume of the WS cell is equal to that of the 
primitive cell, iV ce n/rij. A convenient measure of interparticle spacing is the ion-sphere radius a = (3/47m i ) 1 / 3 . The 
WS cell of the reciprocal lattice is the first Brillouin zone (BZ); its volume equals Vbz = {^Y' n i- 

The frequencies u> s and polarization vectors e s of lattice vibrations (s = 1, . . . ,3Y cc n enumerates vibration modes) 
at any point q of the BZ are determined by (e.g., Ref. ||) 
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■£> a/s (M',q)e?* =0, (1) 

~P fc'/3 

where the summation is over three Cartesian coordinates (Greek indices) and over the ions in the primitive cell (fc'); 

V a0 (k,k', q) = -5 a ,5 kk , - — [g^E (2) 
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is the dynamical matrix, r = l + x(fc) — x(fc'), and x(fc) specifies the ion position within the primitive cell. The primed 
sum means that the term 1 = is excluded if k = k! . 

The elements of the dynamical matrix can be calculated using the Ewald technique of theta-function transformations 
(e.g., Ref. ||), which yields 
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The last sum is over all reciprocal lattice vectors, and p is a free parameter adjusted to yield equally rapid convergence 
of direct and reciprocal sums; a suitable choice is pa w 2. Numerical calculations according to Ecl (||) become unstable 
at qa <C 1 (near the BZ center). In this region we replace D a p by an appropriate asymptote fig ], whose coefficients 
have been recalculated with an accuracy ~ 10~ 8 using the Ewald technique. 
The static-lattice binding energy of a Coulomb lattice is 



E = K M Z 2 e 2 /a, 
where the Madelung constant Km can be written as: 



K M = 
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Previously .Km was calculated, e.g., in Refs. JL7 18|. Our calculated values of .Km for bec, fee, and hep crystals are 
given in Table [j| 



III. BZ INTEGRATION AND FREQUENCY MOMENTS 

In many physical problems, one needs to average functions f(ui) over phonon branches and wave vectors: 

</> = oaT- E /' f=V~ [ W dc *' ( 6 ) 

oiVceU ^ l/_BZ JbZ 

where /(q) = /(w s (q)). In Eq. @ we will use the Holas integration method considered in Ref. [[l9) for the bec lattice: 

/= / d£ /" d»j I dCneHf], (7) 
Jo Jo Jo 

where £, 77, and £ are appropriate BZ coordinates. For the bec crystal, J-{f} — 6/(q), with q = (q x ,q y ,qz) — 
(2 — 77,77,77^)71^/0;, and the lattice constant a; is given by n^af = 2. We calculate the integrals in Eq. (0) by the 

Gauss method involving the nodes of the Jacobi polynomials Pn°'°\ The integral over 77 is alternatively treated by 
the generalized Gauss scheme with weight function 77, which involves the nodes of P^ 1,0 '. 

This approach can be also developed for the fee and hep lattices. In both cases we come again to Eq. (^), but 
with different T{f}. For the fee lattice, we have = §[§/(qi) + §/(q2) + /(<&)], where q 4 = Qi7r£/(2aj), 

Qi = (2 + V + 77C, 2 + V - 77C, 2 - 277), Q 2 = (2 + 277, 2 - 77 + T](, 2 - 77 - T](), Q 3 = (4,ry + t?C, ^ - ??C), and n,af = 4. 

For the hep lattice, T{f} = 2/(q 1 )/jj+2/(q 3 ), where q, = Q, 27r£/(3a,), Qi = (V3, C. f^M Q2 = 7?C, f/ff), 

and riiaf = y/2. Here, a = \/8/3 is twice the ratio of the distance between the hep lattice planes to the distance 
between neighbors within one plane. 

Phonon frequency moments ({lo / Lu p ) n ) and the average (\n(cu / uj p )) , obtained by this method, are given in Table DL 
We remind that ((lo/uj p ) 2 ) = |, according to the Kohn rule (e.g., Ref. |l7]). The accuracy of the data in Table ]l| 
corresponds to the number of digits shown; it is the same or higher than the accuracy of the previous results (e.g., 
Refs. pUp^p^pOj ), except only the value of (lu/lo p ) for the hep lattice, calculated more accurately in Ref. |2C| ] 
((w/cjpVcp = 0.5133368). 
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IV. THERMODYNAMIC FUNCTIONS 



Free energy F of a harmonic Coulomb crystal consists of the static-lattice contribution Eq , contribution from zero- 
point ion vibrations, ^Nh(uj), and thermal free energy in the harmonic lattice approximation, F t ^. Accordingly, the 
reduced free energy / = F j (NksT) is 

f = K M T+ 1.5 (w) + f th , (8) 

where / th (0) = 3 (In (1 - e" 10 )) , and 

aksT ksT ui p 

Thus, the reduced internal energy u = U / {NUbT) = —df/dhxT is 

u = K M r + 1.5 (w) + u th , (9) 

where 

The harmonic constituent of the reduced heat capacity, cy = (Afce) -1 dU/dT = u + du/dlnT, is 

C " (g) ^ Mth -dTn^^ 3 \ (l- e -) 2 /- ( } 

Using the results of Sees. [l] and [II, we have calculated /th(#), u t h(8), and cy(0) for bec and fee crystals as 
corresponding BZ averages. The mean numerical error is estimated as ~ 10~ 6 , and it is a few times larger at 8 3> 1. 
Let us discuss possible analytic approximations. The model of Chabrier et al. [|l3| assumes the linear dispersion law 
for two acoustic (Debye-type) modes, uj± = auj p q/qB 1 and an optical (Einstein-type) mode, u>u — r y-u p . The known 
phonon spectrum moments of a Coulomb crystal are approximately reproduced with the choice a ~ 0.4, 7 ps 0.9. In 
this model, 

/ th = 2 In (1 - e- aB ) + In (1 - e"^) - § D 3 (a8), (12) 

where D%{z) = (3/z 3 ) J^i 3 /^' — l)dt is the Debye function. This model reproduces numerical values of /th, "th, 
and cy with an accuracy of ~ 10%. 

A heuristic generalization of Eq. ( |l2"|) provides a convenient fitting formula to /th- Introducing three logarithmic 
terms (according to three phonon modes) and replacing D3 by an arbitrary rational-polynomial function possessing 
the correct asymptote cx 8~ 3 at large 8, we obtain: 

n=l ^ ' 

where 

8 7 

A(e)=^2a n e n , B(8) = Y / bn0 n + a 6 a 6 6 9 + a 8 a 8 8 1 \ (14) 

n=0 n=Q 

and the parameters a n , a n , and b n are given in Table 0. 

Calculation of the harmonic thermal energy and heat capacity from Eq. (nM using Eqs. @ and @) yields: 



n=l 
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where the first and second derivatives A! \A" ', B' , and B" are readily obtained from Eq. (|14[). 

The approximations (|T^), (15), and (16) have a fractional accuracy within 5 x 1CT 6 , 2 x 1CT 5 , and 5 x 10" 



respectively. 

In the classical limit 9 — ► 0, the exact expansion of fth is 

fa _3 ta . + 3( h (i))-§(£), + ±* + ... (17, 

Note that the term — | (lo/lu p ) 9 cancels the zero-point energy in Eq. (||). Our fit ( |l3| ) reproduces the logarithmic, 
constant, and linear terms of Eq. (|17|) exactly (by construction), whereas the last (quadratic) term is reproduced 
with the relative accuracy of 5 x 10 and 10 -6 for bcc and fee lattices, respectively. Although we do not present 
calculations of the thermal thermodynamic functions for hep crystals, our analysis reveals that they do not deviate 
from the functions for fee crystals by more than a few parts in 10 3 . 

Our results can be used in any applications which require a fast and accurate evaluation of the thermodynamic 
functions of the Coulomb crystals. 
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TABLE I. Parameters of Coulomb crystals. 



lattice type 


K M 


((u,M,)- 2 > 


{(aV^)- 1 } 


((w/wp)) 


<(^M>) 3 > 


(ln(w/Wj,)) 


bcc 


-0.895 929 255 682 


12.972 


2.798 55 


0.5113875 


0.250 31 


-0.831 298 


fee 


-0.895 873 615 195 


12.143 


2.719 82 


0.5131940 


0.249 84 


-0.817 908 


hep 


-0.895 838120 459 


12.015 


2.7026 


0.513 33 


0.24984 


-0.815 97 
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TABLE II. Parameters of the analytic approximation ([L3|) to the thermal free energies of bcc and fee Coulomb lattices. 
Powers of 10 are given in square brackets. 







bcc lattice 






fee lattice 




n 




On 


bn 




a n 


bn 





— 


1 


261.66 


— 


1 


303.20 


1 


0.932446 


0.1839 





0.916707 








2 


0.334547 


0.593586 


7.07997 


0.365284 


0.532535 


7.7255 


3 


0.265764 


5.4814 [-3] 





0.257591 








4 




5.01813 [-4] 


0.0409484 




3.76545 [-4] 


0.0439597 


5 







3.97355 [-4] 







1.14295 [-4] 


6 


4.757014 [-3] 


3.9247 [-7] 


5.11148 [-5] 


4.92387 [-3] 


2.63013 [-7] 


5.63434 [-5] 


7 







2.19749 [-6] 







1.36488 [-6] 


8 


4.7770935 [-3] 


5.8356 [-11] 




4.37506 [-3] 


6.6318 [-11] 
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